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ABSTRACT 

Computer  simulations  were  performed  to  locate  the  equilibrium 
positions  and  binding  energies  of  interstitial  He,  Ne,  Ar ,  Kr ,  and 
Xe  atoms  in  a  tungsten  crystal.   Heavy  interstitial  atoms  in  tung- 
sten share  a  lattice  site  with  the  atom  that  normally  occupies  that 
site  and  form  what  is  called  a  split  interstitial.   Three  charac- 
teristic interstitial  sites  were  located  relative  to  each  lattice 
site  tested.   The  distance  of  the  impurity  atom  from  the  site  was 
seen  to  vary  roughly  inversely  with  its  mass,  and  the  displacement 
of  the  lattice  atom  increased  with  the  mass  of  the  impurity  atom. 

The  foreign  atom  in  its  interstitial  position  was  tested  to 
determine  the  minimum  initial  kinetic  energy  needed  to  escape  the 
lattice,  as  well  as  the  optimum  escape  direction.   The  minimum 
energy  may  be  interpreted  to  be  the  binding  energy  of  the  defect. 
A  comparison  of  experimental  binding  energies  from  Kornelsen  and 
Sinha  and  simulated  binding  energies  indicates  the  model  gives 
realistic  results. 
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I.   INTRODUCTION 

The  advent  of  modern,  high  speed  digital  computers  has  led  to 
the  application  of  computer  simulation  techniques  to  many  differ- 
ent types  of  physical  systems.   One  such  application  is  the  modeling 
of  the  situation  that  occurs  when  a  foreign  atom  interacts  with  a 
metallic  crystal  lattice.   In  general,  such  modeling  can  be  broken 
down  into  two  basic  areas,  dynamic  simulation  and  static  simulation. 
As  an  example  of  the  former ,  radiation  damage  has  been  studied  by 
the  simulated  firing  of  an  atom  or  ion  onto  a  crystal  face.   Other 
examples  are  sputtering  simulations  L1,2,3J  in  which  the  incoming 
particle  causes  surface  atoms  to  be  ejected;  and  channeling  simu- 
lations [4]  in  which  the  ranges  of  ions  travelling  in  crystal 
lattices  are  calculated.   Static  simulations,  on  the  other  hand, 
have  been  concerned  with  the  equilibrium  positions  in  the  lattice 
after  point  defects  such  as  replacement  atoms,  interstitial  atoms, 
and  vacancies  have  been  introduced.   Examples  of  this  type  of 
simulation  can  be  found  in  [5,6,7].   This  present  research  uti- 
lized aspects  of  both  static  and  dynamic  simulation  techniques. 
The  goal  of  this  research  was  to  correlate  the  results  of  experi- 
mentally determined  binding  energies  of  point  defects  in  a  tungsten 
lattice  [8,9]  with  the  results  obtained  by  computer  simulation. 


II.   NATURE  OF  THE  PROBLEM 

A.   HISTORICAL  BACKGROUND 

An  investigation  of  radiation  damage  events  by  computer  simu- 
lation techniques  of  crystalline  behavior  was  published  in  i960  by 
Gibson,  Goland,  Milgram,  and  Vineyard,  hereafter  referred  to  as 
(GGMV)  [5].   This  Brookhaven  National  Laboratory  investigation  set 
forth  the  criteria  that  must  be  satisfied  before  the  simulation 
method  could  be  applied  to  crystals.   Such  factors  as  potential 
energy  functions,  forces,  crystallite  sizes,  computation  methods, 
choice  of  time  intervals,  and  computer  limitations  were  discussed. 
The  crystal  lattice  modeled  in  their  research  was  metallic  copper, 
a  face-center  cubic  (fee)  structure.   The  potential  functions  used 
in  the  calculations  was  a  Born-Mayer  repulsive  potential,  with  the 
necessary  cohesive  forces  applied  on  the  boundries  of  the  crystal- 
lite.  In  integrating  the  equations  of  motion,  the  Brookhaven 
group  used  the  central  difference  procedure.   The  optimum  choice 
for  timestep  duration,  At,  was  shown  to  be  of  critical  importance 
in  the  integration  scheme.   The  energy  of  the  strongest  interac- 
tion governed  their  choice  of  the  above  parameter.   The  static 
results  obtained  by  GGMV  confirmed  the  existance  of  the  (l00> 
split  interstitial  in  the  fee  lattice.   Their  dynamic  results 
described  collision  chains  and  focusing  phenomena  in  crystallites 
struck  with  energetic  knock-on  atoms. 

Additional  crystal  simulation  studies  were  performed  by 
R.A.  Johnson  L6,10,11J.   In  Ref.  6,  he  investigated  point  defects 
in  a  copper  lattice  using  Born-Mayer  repulsive  potentials.   The 


^10(3)  split  interstitial  was  found  to  be  the  only  stable  inter- 
stitial position.   He  found  it  necessary  to  allow  the  interstitial 
to  interact  only  with  its  six  nearest  neighbors.   Atoms  near  the 
defect  were  treated  as  independent,  while  the  remainder  of  the 
metal  was  treated  as  an  elastic  continuum  with  atoms  imbedded  in 
it. 

Research  on  body-centered  cubic  (bcc)  crystals  was  undertaken 
by  Erginsoy,  Vineyard,  and  Englert  (EVE)  L7l.   They  used  a  com- 
posite potential  function  for  most  of  their  detailed  work.   It 
consisted  of  an  exponentially  screened  Coulomb  potential  at  small 
separations,  a  Born-Mayer  function  in  the  region  between  small 
and  intermediate  separations,  and  a  Morse  potential  at  larger 
separation  distances.   A  split  interstitial  was  reported  for  simu- 
lated bcc  crystals  in  the  (lio)  direction.   In  their  dynamic  re- 
sults they  reported  the  existence  of  a  threshold  energy  for 
displacement  that  was  highly  independent  of  the  direction  of 
knock-on.   Also,  collisional  chains  in  the  (ill)  and  (lOO/  direc- 
tions were  found  to  exist.   R.A.  Johnson  also  published  results 
for  bcc  simulations  involving  a-iron  and  tungsten  LllJ.   The 
existence  of  the  split  interstitial  in  the  (lio)  direction  was 
confirmed  and  crowdian  migration  data  were  discussed.   The  present 
investigation  confirmed  the  (lio)  split  interstitial  positions 
for  argon,  neon,  krypton,  and  xenon. 

D.E.  Harrison  and  associates  have  published  several  articles 
in  which  computer  simulation  of  crystalline  behavior  has  been 
investigated  [l,2,12].   In  a  study  of  a  fee  model  of  copper,  col- 
lision events  between  a  copper  atom  and  a  copper  lattice  were 
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simulated  as  a  function  of  the  potential  function,  the  enegy  of 
the  collision,  and  the  location  of  the  impact  point.   The  inte- 
gration scheme  used  was  an  average  force  procedure,  instead  of 
the  central  difference  procedure  used  by  GGMV.   A  complete  dis- 
cussion of  the  method  can  be  found  in  Ref.  13.   A  continuation  of 
copper  simulations  was  published  by  Harrison,  Leeds,  and  Gay  in 
1965  [12].   Another  paper  by  Harrison  and  Greiling  dealt  with  the 
ranges  of  heavy  ions  in  tungsten  crystals  whose  atoms  had  under- 
gone thermal  displacement  L4J-   It  was  found  that  room  temperature 
thermal  displacements  had  a  negligible  effect  upon  the  collisions 
for  ion  energies  greater  than  a  few  thousand  electron  volts. 
Finally,  Harrison,  Levy,  Johnson,  and  Ef f r on  published  results  on 
computer  simulation  of  sputtering  [2J . 

Research  undertaken  by  Vine  Ll4J  provided  the  foundation  for 
this  author's  present  investigation.   Vine  used  the  Gay-Harrison 
model  for  crystal  simulation  as  modified  by  Levy  Ll5],  Johnson  Ll6J, 
Effron  [17],  and  Moore  Cl8] .   His  overall  objective  was  to  cor- 
relate experimental  and  simulated  binding  energies  of  neon  and 
argon  point  defects.   Repulsive  potential  functions  for  neon- 
tungsten  (Ne-W)  and  argon-tungsten  (Ar -W)  interactions  were 
used  [19].   Morse  functions  were  not  used  for  those  interactions 
because  experimental  data  giving  the  Morse  parameters  was  based  on 
homogeneous  media  such  as  tungsten-tungsten  (W-W)  interactions  L20J. 
Vine  subsequently  attempted  to  correlate  results  of  equilibrium 
position  studies  for  tungsten  defects  in  a  tungsten  lattice  using 
two  different  potential  functions  for  the  interstitial-lattice 
interactions.   In  one  case,  the  tungsten  interstitial  was  allowed 
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to   interact    using   a    Born-Mayer    repulsive  potential,    and   the   other 
case,    the    tungsten    defect   was    given   a    composite  potential    that   was 
identical    to   that   given    to  all    the    other    lattice   atoms.       By    so 
doing,    he   attempted    to   establish   an    empirical    relationship   between 
the   two  methods    to   apply    to   the   results    of   neon    interstitial    studies 
which   used    only   a    repulsive  potential. 

The   results    of   the    tungsten-tungsten    interactions    failed   to 
provide   the   information    needed   to   formulate   a    correction    factor 
to   be   used   in    the   neon-tungsten    studies.       In   addition,    the   concept 
of   relating    the  potential    energy   at    equilibrium    to   the    experi- 
mentally   observed   binding    energies    of   Kornelsen    and  Sinha    [9]    does 
not   appear    to   be   feasible. 

B.       CHOICE   OF   POTENTIAL    FUNCTION 

The    studies    that   were   reviewed   in    the  previous    section   utilized 
many    different    approximations    to   the   true  potential    function    be- 
tween  atoms    in   a    metal    lattice.      The  problem   of   solving    the   many- 
body    interaction    of  a    real    system    is   approximated    in    the   computer 
by   many    two-body    interactions.      Thus,    central   pairwise  potential 
functions    are   most    often   used    in   computer    simulations.      GGMV 
employed  a   repulsive  potential    of   the   Born-Mayer    form: 

V.  .    =    exp(A+Br .  .) 

which   described   the   repulsion    of  atoms    at   close   approach.      Three 
Born-Mayer    potentials    were    investigated  by   GGMV   to    determine   which 
would   give   the  best   results    in   their    calculations.      Their    choice 
was    one   which   has    since   been    labeled   the  Gibson   Number    Two  Potential 
In   Ref.    6,    Johnson   and   Brown    used  a    similar    Born-Mayer    potential 
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with  slightly  different  parameters.   As  previously  mentioned, 
another  potential  function  that  has  been  used  in  crystal  modeling 
studies  is  the  familiar  Morse  potential  of  the  form: 

*ij  =  D[e*P  ["  2a(rij  "  re1j  "2  «*{■*&«  "  re)}] 

where   r      is    the   equilibrium    distance    of   approach    of   two  atoms, 
and  Ci   and   D  are    constants. 

Girifalco  and  Weizer    calculated  Morse   parameters    that    would 
be  appropriate   for    several    crystal    lattices    [20].      In    calculating 
the  parameters,    they   attempted   to   express    the    various    crystal 
properties    such   as    cohesive   energy,    lattice    constant,    compres- 
sibility  and  equation    of   state    in    terms    of   the  Morse    function. 
The  Morse   potential    constants   published  by   Girifalco   and  Weizer 
have   been    used  extensively    in    simulating    the  potential    functions 
and  forces    in    lattices    of   homonuclear    atom   systems. 

The    Born-Mayer    potential   and   the  Morse  potential   are   useful 
over    specific    internuclear    separation    distances.      The    Born-Mayer 
potential    is    useful   at    strongly   repulsive   separations,    i.e.    short 
ranges,    while    the  Morse    function    is    applicable   at    equilibrium  and 
greater    separations.       In    order    to   better    approximate    the    true 
potential    function,    various    investigators    have    combined    different 
potential    functions.      As   mentioned   earlier,    EVE   [7]    used  a    com- 
posite potential    that    consisted   of   an    exponentially    screened 
Coulomb  potential   at    very    close   separations,    a    Born-Mayer    po- 
tential  at   weakly   repulsive   distances,    and  a    Morse  potential    over 
the   remainder    of   the  potential    curve.       In    their    studies,    Harrison 
and  associates    combined  a    Born-Mayer    potential   and  a    Morse 
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potential.      These   two  potentials    were   joined  by  a   cubic    equation 
in   the  region   near    their    intersection   Ll,2,3,4l. 

C.       EXPERIMENTAL    RESULTS   OF  KORNELSEN   AND  SINHA 

In   1968,    Kornelsen   and   Sinha    L7,8]    published  results    concerning 
the   binding   energies    of   trapped  particles    such   as    neon,    argon, 
krypton,    and  xenon    ions    in   a    tungsten    surface.      The  particles    were 
forced   into   the   surface   from  a   beam   created   by   an    ion   gun  which 
gave    ion   energies    of  40    eV   to   5   keV.      The   tungsten   crystal   was 
then   heated  and   the  rates    of   evolution   of   the   trapped  gas    were 
measur  ed. 

The   temperature   at   which    desorption   peaks    occured   thus    gave 
an    indication   of   the   binding    energies    of  point    defects    in    the 
tungsten    crystal.       Quantitatively    siroiliar    results    were    obtained 
with   argon,    krypton,    neon,    and   xenon    ions.      Four    peaks    were    ob- 
served below    1650      K   and  were    labeled  as   0!  peaks.      A   single   peak 
above    1700      K   was    measured  and   was    called   the   p-peak.       It   was    con- 
cluded  that    the  a-group    of  peaks    were   the   results    of   single   step 
desorptions    from   sites   very    close   to    the   metal    surface.      They 
further    concluded   that    the   different  tt-peaks    could   correspond   to 
different    types    of  point    defect    binding    energies    in    the   tungsten 
crystal.      Specifically   mentioned  were    defects    of    three   types; 
(a)    interstitial   and   substitutional   positions    in    the    lattice,    (b) 
different    distances    from   the    site   to  the   surface,    and    (c)    different 
locations    of  a   nearby    lattice    defect,    such   as   a    vacancy,    relative 
to  the   site  and  the   surface. 
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III.       THE   SIMULATION   MODEL 

A.       THE   CRYSTAL 

The  model  used  in  this  research  was  essentially  the  same  one 
used  by  Vine  Ll4J  as  explained  in  Section  II.   In  subsequent  dis- 
cussion, when  it  is  necessary  to  specify  certain  of  the  computer 
program  variable  names,  they  will  be  placed  in  braces. 

All  of  the  simulations  in  this  research  were  done  with  tungsten 
crystals  of  varying  sizes.   The  objective  was  to  use  the  smallest 
crystal  dimensions  that  would  give  realistic  results.   The  di- 
mensions described  below  refer  to  the  number  of  planes  of  atoms 
in  each  of  the  three  rectangular  coordinate  directions.   Simu- 
lations were  performed  with  sizes  8  x  6  x  8  of  96  atoms, 
10  x  6  x  10  of  150  atoms,  10  x  8  x  10  of  200  atoms,  and 
10  x  10  x  10  of  250  atoms.   Of  these,  the  latter  two  were  judged 
to  be  of  most  use  because  of  their  greater  depth  in  the  y-direction. 
The  y  direction  was  always  used  as  the  direction  of  escape  for  the 
point  defect  atom. 

Each  atom  in  the  simulated  crystal  was  numbered,  with  the 
first  position  always  assigned  to  the  point  defect  atom.   The 
remainder  of  the  positions  were  assigned  in  sequence  according  to 
the  coordinate  locations.   The  numbering  was  started  in  the  y  =  0 
plane  and  continued  until  all  the  atoms  in  that  plane  were  speci- 
fied.  This  procedure  was  repeated  for  the  remainder  of  the  y  planes 
in  the  simulated  crystal. 
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B.   THE  TIMESTEP  INTERVAL 

The  numerical  method  of  time  integration  used  in  the  model  was 
the  average  force  method  Cl3l.   The  value  of  At  {dt}  used  in  this 
procedure  was  of  critical  importance  in  determining  whether  or  not 
the  model  would  approximate  reality.   Also,  the  program  running 
time  was  a  function  of  the  timestep  duration. 

In  order  to  best  approximate  reality,  the  iDTJ  was  kept  smaller 
early  in  the  program  when  most  movement  was  expected,  and  was  al- 
lowed to  grow  larger  as  equilibrium  was  approached.   The  parameter 
that  controlled  the  timestep  duration  was  {dti},  the  distance  which 
the  most  energetic  atom  was  allowed  to  move  in  a  single  timestep. 
In  previous  work,  {dti}  was  held  constant  throughout  the  duration 
of  a  program.   If  the  motion  was  expected  to  be  slow,  IDTI]  was 
given  a  larger  value  than  in  a  situation  where  simulated  motion 
was  expected  to  be  greater.   For  static  equilibrium  problems,  at 
least  100  timesteps  were  needed  to  reach  an  approximately  stable 
position.   In  addition,  it  became   necessary  to  insert  a  maximum 
value  for  the  timestep,  { Dt} ,  into  the  program  to  prevent  unrea- 
listic movement  of  the  atoms  and  breakdown  of  the  model. 

In  the  current  program,  changes  were  made  to  allow  {dTIj  to 
decrease  during  the  program.   For  example,  in  static  runs  on  He-W, 
Ne-W,  Ar-W,  and  Kr-W,  [dti}  was  initially  given  a  value  of  0.1 
lattice  unit  {lu} .   Tungsten  forms  a  bcc  crystal,  with  LU  equal 
to  ?gLC  or  1.58  A,  where  LC  is  the  Lattice  Constant  or  cube  edge 
distance.   {CTl}  was  allowed  to  change  in  decrements  of  0.01  LU 
per  timestep  for  the  first  10  timesteps.   Then  the  {dti}  value 
of  0.01  LU  was  allowed  to  change  in  decrements  of  .001  LU  for 
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another  10  timesteps.   Equilibrium  was  reached  by  30  timesteps 
according  to  this  procedure,  resulting  in  a  significant  decrease 
of  computer  running  time.   Also,  the  equilibrium  positions  obtained 
in  the  simulation  were  closer  to  the  expected  (lio)  split  inter- 
stitial positions  than  were  obtained  with  a  constant  iDTl}  value. 
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IV .   SIMULATION  PROCEDURE 

A.   STATIC  SIMULATIONS 

The  first  stage  of  the  simulation  procedure  was  concerned  with 
finding  the  equilibrium  positions  for  the  point  defect  of  interest 
in  the  top  several  layers  of  a  bcc  tungsten  crystal.   Previous  work 
done  by  Johnson  [ll],  indicated  that  a  bcc  crystal  had  at  least 
two  different  types  of  stable  split  interstitial  orientations. 
The  first  of  these  was  in  a  (lOO)  direction  along  a  (110)  plane. 
The  other  stable  orientation  was  in  a  (lio)  direction  in  a  (100) 
plane.   An  analysis  of  the  total  of  12  stable  split  interstitial 
positions  possible  about  a  given  atom  in  the  two  orientations 
listed  above,  revealed  that  there  were  only  three  independent  types 
of  sites.   (See  Figure  1.)   The  first  of  these  is  located  on  a 
(110)  plane  with  {nVAC}  and  lies  closer  to  the  surface  than  {nVAc} . 
The  second  independent  position  lies  in  the  same  (100)  plane  as 
iNVACj  and  is  at  the  same  depth  in  the  crystal.   A  final  position 
similar  to  the  first,  in  that  it  is  along  a  (110),  is  located  at 
a  depth  below  that  of  {nVAC}.   The  three  different  split  inter  - 
stitials  were  postulated  to  have  different  binding  energies  be- 
cause of  their  varying  depths  in  the  crystal. 

In  order  to  reduce  the  running  time  in  locating  the  final 
coordinate  positions  at  each  interstitial  location,  the  static 
programs  were  started  with  the  foreign  atoms  in  the  approximate 
area  of  the  expected  final  positions  as,  discussed  above.   For 
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Figure   1 

SPLIT    INTERSTITIAL       SITES     FOR  BCC     TUNGSTEN 
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example,  when  it  was  desired  to  locate  position  A,  the  foreign 
atom  was  initially  placed  along  the  correct  (110)  plane  approxi- 
mately one  lattice  unit  from  {nVACj. 

B.   DYNAMIC  SIMULATIONS 

The  dynamic  simulation  program  used  the  final  positions  com- 
puted in  the  static  programs  as  input  initial  crystal  positions. 
The  lattice  generator  subroutine  { BIOO} ,  and  the  point  defect 
locator  {PLACE}  were  thus  eliminated  from  the  dynamic  program. 

Dynamic  simulations  were  broken  down  into  two  main  categories. 
The  first  category  consisted  of  survey  runs  of  the  areas  above 
the  interst itials  in  the  direction  of  escape  from  the  y  surface. 
An  impact  point  generator  package  was  inserted  into  the  dynamic 
program,  thus  permitting  the  interstitial  to  be  directed  at  a 
specific  number  of  locations  in  a  predefined  area.   The  results 
of  the  survey  run  were  analyzed  to  determine  the  optimum  aiming 
point  for  the  interstitial  at  a  specific  initial  energy.   It  was 
observed  during  the  impact  testing  precedure  that  the  "best" 
aiming  points  were  a  function  of  the  initial  energy  given  to  the 
interstitial.   However,  it  was  also  found  that  the  optimum  aiming 
points  at  varying  energies  were  in  a  generally  localized  area. 
Thus,  once  the  optimum  points  were  found  at  the  starting  energy, 
only  a  localized  region  around  those  points  was  tested  at  lower 
initial  interstitial  energies.   The  decrementing  process  of  the 
initial  interstitial  energy  constituted  the  second  main  phase  of 
the  dynamic  simulations.   The  procedure  in  this  phase  was  to  de- 
crease the  interstitial  energy  until  it  could  no  longer  escape 
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from  the  crystal.   The  minimum  escape  energy  was  said  to  be  the 
simulated  binding  energy  of  the  interstitial  at  the  particular 
location  tested. 
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V.   PRESENTATION  OF  CATA 

A.   STATIC  RUNS 

1 .   Preliminary  Testing  of  Static  Program 

One  hundred  twelve  computer  runs  were  made  in  the  initial 
testing  phase  of  the  static  program.   The  bulk  of  the  testing  was 
done  in  four  general  areas:  1)  Program  shutdown  procedures  at 
equilibrium;  2)  Optimum  crystal  size  determination;  3)  Equilibrium 
position  as  a  function  of  the  initial  interstitial  position;  and 
4)  Realistic  timestep  determination  procedure. 

The  best  method  of  stopping  the  program  at  equilibrium  has 
been  a  matter  of  concern  for  many  years  with  the  static  simulation 
program.   In  the  present  investigation,  the  first  method  tested  was 
a  shutdown  procedure  initiated  whenever  a  sharp  { DT}  decrease  was 
encountered  in  the  program.   This  test  proved  to  be  of  limited 
success,  and  was  later  abandoned.   Another  method  that  was  at- 
tempted was  a  test  of  [eMAX]  against  a  value  such  as  .04,  to  deter- 
mine if  equilibrium  had  been  reached.   This  method  also  proved  only 
partially  successful.   Finally,  a  test  was  made  of  the  average 
kinetic  energy  of  an  atom  in  the  simulated  crystal.   The  average 
energy  was  taken  to  be  the  total  kinetic  energy  {tPKe}  divided  by 
the  number  of  atoms  {ll} ,  or  equivalently ,  {TPKe}  multiplied  by 
the  reciprocal  of  the  number  of  atoms  { RLLJ .   The  crystal  was 
assumed  to  be  at  equilibrium  if  the  average  kinetic  energy  of  an 
atom  was  less  than  or  equal  to  the  value  0.025  eV,  a  value  for  the 
average  thermal  energy  of  an  atom.   Satisfactory  results  were 
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sometimes  obtained  by  this  method  but  the  test  had  to  be  removed 
in  many  cases  to  allow  the  simulated  crystal  to  run  a  greater 
number  of  timesteps  and  reach  equilibrium. 

As  previously  mentioned,  different  sizes  of  crystals  were 
simulated  to  determine  the  optimum  dimensions  for  the  crystallite. 
Many  tests  were  made  on  crystal  sizes  smaller  than  the  10  x  10  x  10 
used  by  Vine  Cl4] .   Although  smaller  sizes  such  as  10  x  6  x  10  often 
gave  reliable  results,  it  was  finally  decided  to  use  the  10  x  10  x  10 
size  for  the  reported  split  interstitial  positions  in  order  to  have 
a  standard  size  applicable  over  all  crystal  positions  of  interest. 

In  his  work,  Vine  placed  his  inter stitials  in  the  middle 
of  the  open  channels  in  the  simulated  crystal.   Numerous  runs  in 
the  present  research  have  indicated  that  equilibrium  is  reached 
more  quickly  when  the  initial  starting  positions  are  not  in 
channel  centers,  but  along  the  directions  of  the  expected  splits 
in  the  general  area  of  the  final  positions. 

All  of  the  initial  work  was  done  with  a  fixed  {CTl}  value 
in  the  program.   The  [dti]  value  normally  used  was  .05  LU.   The 
timestep  interval  was  later  modified  as  explained  in  The  Simulation 
Model  and  the  computer  running  time  was  cut  by  approximately  two- 
thirds  due  to  this  procedure. 

2.   Positions  for  Helium  in  Tungsten 

Helium  was  the  lightest  point  defect  used  in  the  static 
simulation  model.   It  was  thought  that  the  small,  light  atom  would 
essentially  do  all  the  movement  and  come  to  rest  in  a  channel 
center  \J~2~  LU  away  from  [iWAc}  along  the  (ll0>  direction.   The 
results  of  the  runs  show  that  the  movement  was  almost  as  expected. 
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Also,  a  comparison  of  the  C  site  for  atom  64  and  the  A  site  for 
atom  114  indicates  that  there  is  a  strong  possibility  that  these 
two  sites  are  degenerate.   This  is  based  on  the  fact  that  the 
potential  energies  of  both  sites  are  close,  and  also  that  the 
{nVAc}  for  each  site  is  displaced  only  a  negligible  distance.   The 
same  probability  of  degeneracy  is  also  seen  to  exist  for  C-89  and 
A-139.   The  determination  of  degeneracy  of  sites  is  seen  to  be  a 
complex  evaluation  of  potential  energy  differences,  movements  of 
{nVAC} ,  and  the  potential  gradient  between  the  two  sites.   The 
final  positions  for  the  A,  B,  and  C  split  interstitial  positions 
in  the  third,  fourth,  fifth,  and  sixth  layers  where  {NVAC}  was 
64,  89,  114,  and  139  respectively  are  shown  in  Figures  2,  3  and  4. 
The  results  obtained  with  helium  were  thought  to  be  somewhat  in 
error,  since  the  heavier  atom,  neon,  moved  further  in  its  simu- 
lations than  helium.   Further  work  is  needed  to  confirm  the  equi- 
librium positions  for  helium. 

3.  Positions  for  Neon  in  Tungsten 

The  neon  split  interstitial  locations  were  simulated  for 
the  same  positions  as  helium.   (See  Figures  5,  6  and  7.)   Essenti- 
ally, the  [NVAC]  atom  remained  in  its  lattice  site  and  the  neon 
moved  along  a  (lio)  direction  to  a  distance  \|2  LU  away  from 
1  NVAC}  . 

4.  Positions  for  Argon  in  Tungsten 

The  bulk  of  the  testing  of  this  research  was  done  with 
argon  as  the  simulated  point  defect.   The  results  of  the  runs  for 
the  10  x  10  x  10  crystal  size  are  shown' in  Figures  8,  9  arid  10. 
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Figure   3 
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Figure   4  
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Figure  6 

A    a    C    SITES   FOR    ATOMS     89,      l~3~9 

TUNGSTEN  -  O 


2  n*»  NEON   - 


3.0 


(Y)LU    4.0 


\ 


\ 


\ 


\ 


\ 


/ 


/ 


/ 


/ 


/ 


s*/ 


5,0 


\ 


\ 


\ 


\ 


\ 


\ 


\ 


\ 


\ 


\ 


\ 


\ 


/ 


/ 


/ 


/ 


/ 


/ 


/ 


/ 


139} 


/ 


/ 


/ 


/ 


/ 


6.0 


\ 


\ 


\ 


\ 


\ 


\ 


\ 


4.0 


5.0 
(X)    LU 


6.0 


29 


Figure   7 

B      SITES    FOR    ATOMS    64,89,   114,13a 
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Figure    8 
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Figure   9 
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Figure    10 

B     SITES     FOR    ATOMS    64,    89,114,139 
ARGON  -@  TUNGSTEN  -Q) 
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The  expected  (lio)  splits  were  observed  with  [nvACj  moving  away 
from  its  site  in  relation  to  the  final  position  of  the  argon  defect 
5.   Positions  for  Krypton  in  Tungsten 

The  split  interstitial  positions  for  krypton  in  a  simu- 
lated tungsten  crystal  were  also  plotted.   (See  Figures  11,  12 
and  130   Once  again,  satisfactory  results  were  obtained  and  the 
(llO/  splitting  was  observed. 

^  *   Positions  for  Xenon  in  Tungsten 

For  the  Xenon  Runs,  the  tungsten  repulsive  potential 
function  was  used  for  the  xenon  potential.   The  initial  runs, 
using  the  iDTl}  employed  for  the  other  defects,  failed  to  give 
the  postulated  split  interstitial  positions.   At  least  two  factors 
were  seen  to  complicate  the  Xenon-Tungsten  runs.   Firstly,  the 
mass  of  Xenon  was  approximately  6/7  that  of  tungsten,  so  the 
lattice  site  would  be  shared  almost  evenly.   This  would  mean  that 
tungsten  would  have  to  move  a  significant  distance  from  its  lat- 
tice position.   Secondly,  the  relatively  large  size  of  the  xenon 
defect  would  require  more  movement  of  the  surrounding  atoms  to 
accommodate  the  extra  atom. 

The  {DTI}  scheme  was  changed  to  allow  for  more  movement  of 
the  atoms  by  allowing  iDTl}  to  change  in  smaller  decrements. 
Another  method  that  was  used  was  an  initial  displacement  of  both 
the  xenon  and  the  tungsten  away  from  the  lattice  site.   When  both 
atoms  were  displaced,  it  became  necessary  to  use  an  initial  [DTl} 
of  .05  LU  or  less.   An  average  value  of  the  C  site  for  atom  89 
was  computed.   (See  Figure  14 . )   More  work  is  needed  to  simulate 
all  of  the  interstitial  positions. 
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Figure    12 
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Figure    13 

B      SITES     FOR    ATOMS      64,   89,  114,  139 
KRYPTON  —Q  TUNGSTEN  -  O 
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Figure    14 
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7 .   Split  Interstitial  Distance  Ratios  as  a  Function  of 

Relative  Masses 

During  the  static  testing  phase,  the  ratios  of  the  split 
interstitial  distances  from  the  initial  [nvACj  were  measured.   An 
attempt  was  made  to  correlate  these  ratios  to  the  inverse  ratios 
of  the  atomic  masses  of  the  atoms  involved.   (See  Table  I.)   The 
point  defect  results  of  helium  and  neon  did  not  give  any  signi- 
ficant correlation,  but  the  Argon  and  krypton  atoms  gave  ratios 
of  splits  in  good  agreement  with  the  expected  values  from  mass 
ratio  calculations. 

B.   DYNAMIC  RUNS 

1 .   Survey  of  the  Possible  Directions  of  Escape  for  the  Ar -W 
Simulated  Split  Interstitial  in  Site  A-89 

As  was  mentioned  in  the  Section  IV,  survey  runs  were  made 
of  the  possible  escape  directions  in  the  plane  one  LU.  above  the 
defect.   The  results  for  the  survey  for  the  Ar -W  split  interstitial 
in  site  A-89  were  plotted.   (See  Figure  15-)   I DYj ,  the  distance  in 
lattice  units  that  the  argon  moved  in  25  timesteps,  and  {  VYj  ,  the 
velocity  that  the  argon  had  after  25  timesteps  are  shown  for  each 
impact  point.   Due  to  the  multiple  scattering  of  the  defect  as  it 
moved  toward  the  surface,  25  timesteps  were  not  sufficient  to 
provide  conclusive  evidence  at  any  one  impact  point.   The  survey 
was  limited  to  25  timesteps  per  impact  point,  however,  due  to 
computer  time  considerations.  The  main  benefit  gained  from  the 
survey  run  was  the  elimination  of  certain  areas  from  further  testing, 
such  as  those  on  the  outer  perimeter  of  the  survey  area.   Also,  much 
information  was  gained  on  the  most  likely  area  of  escape  for  the 
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TABLE    I 


Theoretical    Splits 
From   Lattice   Site 

^Tungsten   Mass        ), 
\   Impurity  Mass      / 

Simulated  Splits     From 
Lattice   Site 

/    Impurity    Distance\ 
I  Tungsten    Distance      J 

He-W 

45.96 

Ne-W 

9.11 

Ar-W 

4.60 

4.36 

Kr-W 

2.19 

2.18 

Xe-W 

1.40 
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Figure   15 

IMPACT  AREA:  CMIAMIC   SIMULATION  OF  ESCAPE 
DIRECTIONS    FOR    ARGON    IN     SITE      A-  89 
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defect.   From  the  results  of  the  survey  run,  and  from  symmetry 
considerations  of  the  open  channel,  it  was  decided  to  concentrate 
further  testing  on  several  points  with  the  same  Z  coordinate  as 
the  interstitial  atom. 

2 .   Determination  of  the  Simulated  Binding  Energy  of  Argon 

in  Site  A-89 

Detailed  testing  of  several  impact  points  with  the  same  Z 
coordinate  as  the  argon  was  carried  out  for  100  timesteps  per 
impact  point.   This  was  generally  enough  time  for'  the  simulated 
interstitial  to  escape  the  crystal  it  its  initial  energy  was 
great  enough.   The  results  are  summarized  below. 

TABLE  II     Detailed  Impact  Point  Testing 


Impact  Points 
Plane  1.0  LU 
-y  direction 
Argon 

in 
in 

from 

; 

[nitial  Energy 

20  eV 

10  eV 

8  eV 

6  eV 

4  eV 

cx  =  3.84  LU 
cz  =  4.99 

Yes 

Yes 

No 

No 

No 

cx  =  4.24 

CZ  =  4.99 

Yes 

Yes 

Yes 

Yes 

No 

cx  =  4.44 
CZ  =  4.99 

Yes 

Yes 

Yes 

Yes 

No 

Yes  -  Atom  Escapes  Crystal 
No   -  Atom  Does  Not  Escape 
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VI.   CONCLUSIONS 

The  method  used  to  decrease  [DTI]  on  every  timestep  appears 
to  be  a  significant  improvement  over  the  old  method  of  keeping 
{DTI}  fixed.   In  effect,  the  new  method  forces  the  simulated 
crystallite  to  come  to  equilibrium  in  a  predefined  number  of 
timesteps.   In  addition  to  a  saving  of  computer  time,  this  method 
also  eliminates  the  need  for  a  equilibrium  shut  down  procedure 
in  the  program.   The  same  iDTl}  decrement  scheme  could  not  be  used 
for  all  the  point  defects  tested  in  this  research.   For  example, 
Argon  and  Krypton  were  able  to  come  to  their  expected  equilibrium 
positions  using  the  scheme  described  in  Section  III,  but  Xenon 
and  Helium  were  not.   Thus  it  appears  that  defect  size,  as  well 
as  the  expected  degree  of  movement  may  dictate  the  [dTIJ  decre- 
ment scheme. 

The  expected  (lio)  split  inter stitials  were  confirmed  for  the 
helium,  neon,  argon,  krypton,  and  xenon  foreign  defects.   The  re- 
lative degree  of  splitting  does  appear  to  be  related  to  the  masses 
of  the  interacting  atoms,  but  conclusive  evidence  of  this  was  only 
obtained  for  the  Argon-Tungsten,  and  Krypton-Tungsten  runs. 

The  dynamic  runs  have  shown  the  direction  of  escape  from  the 
crystal  to  be  a  mult i-collis ional  process,  with  the  defect  under- 
going many  intermediate  direction  changes  at  the  low  kinetic 
energies  tested.   The  most  likely  escape  direction  was  also  seen 
to  vary  somewhat  with  the  initial  kinetic  energy  given  to  the 
defect . 
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The  simulated  value  of  approximately  4  eV  for  the  binding 
energy  of  Argon  in  site  A-89  appears  to  be  a  reasonable  value. 
Kornelsen's  data  showed  an  energy  level  at  this  value  and  his 
levels  were  postulated  to  arise  from  defects  in  the  first  several 
layers  of  the  tungsten  crystal. 

Further  work  is  needed  to  confirm  the  remainder  of  the  energy 
levels  found  by  Kornelsen,  by  testing  other  point  defect  sites  in 
the  top  layers  of  the  simulated  crystal. 
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APPENDIX  A:  COMPUTER  PROGRAM  GLOSSARY 


NOTE:   In  this  glossary,  the  terms  "point  defect  atom",  "bullet", 

and  "primary"  are  synonymous;  and  the  terms  "latttice  atom"  and 

"target"  are  synonymous. 

AC:       Distance  measurement  used  in  impact  point  generator 

ALPHA:    Input  Morse  potential  parameter 

BSAVE:    Target  mass/(target  mass  +  bullet  mass);  distributed 

potential  energy  between  target  and  bullet 
BIND:     Negative  of  the  total  potential  energy  (TPOT)  at  time 

zero 
BMAS:     Mass  of  bullet  in  amu 

BULLET:   Alpha -numer ic  array  for  point  defect  material 
CFO,  CFl ,  CF2 :   Force  parameters  of  cubic  fit  between  Morse  and 

Born-Mayer  functions 
CGB1 ,  CGB2 :   Morse  potential  parameters 
CGD1 ,  CGD2  :   Morse  potential  parameters 
CGF1 ,  CGF2 :   Morse  force  parameters 

COX,  COY,  COZ:   Cosines  of  angles  to  x,y,  and  z  axes  respectively 
CPO,  CP1 ,  CP2 ,  CP3:  Potential  parameters  of  cubic  fit  between 

Morse  and  Born-Mayer  functions 


CVD:  CVR  x  10 


-10 


converts  lattice  units  to  meters 


,-19 


CVE:    1.6   x    10         ,    converts    electron    volts    to   joules 

CVED:    CVE/CVD,    a    ratio  used   to   avoid  repeated   division 

-27 
CVM:    1.672   x   10         ,    converts    atomic   mass    units    to   kilograms 

CVR:    LU    in   angstroms;    converts    lattice    units    to  angstrom  units 
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CX,    CY,    CZ:       Coordinates    of    impact   point 

D1X,    D1Y,    DlZ:      Displacement    coordinates    for    location    of    inter- 
stitial   from  reference  atom,    NVAC 
DCON:  Input   Morse  potential   parameter 

DDTI:  Time    increment    that    is    subtracted  from    DTI   after    each 

timestep 
DFF:  ROE-DIST,    the    distance   closer    than   ROE   that   an   atom    is    to 

the  pr  imar  y 
DIST:  Distance   between   any    two  atoms 

DLPE:  TLPE-TLPE0,    the    change    in    total   local   potential    energy 

since   time  zero 
DRX,    DRY,    DRZ:    x,y,z    components    of    DIST 
DT:  Length   of  a    timestep    in   seconds 

DTI:  Number    of    lattice   units   most    energetic   atom  may   move    in 

one    timestep 
DTIS:  Starting    value    of    DTI. 

DTOD:  DT/CVD--a    ratio  used    to  avoid   repeated   division 

DTOM:  DT/PTMAS--a    ratio  used   to  avoid  repeated   division 

DTOMB:         DT/PEMAS--a   ratio   used  to  avoid  repeated   division 
DX(I),    DY(I),    DZ(I):      Change    in   position    of    ith   atom   from    initial 

position   at    time   zero 
EMAX:  The   maximum  energy    encountered    in   any    cycle 

ERAT:  Measure   of   the   average    kinetic    energy   of  an   atom 

EV:  Primary    energy    in    electron    volts 

EVR :  Primary    energy    in    kilo-electron    volts 

EXA,    EXB:       Input    Born-Mayer    potential    function   parameters    for    the 

target 
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F2 :  Square   of    the    force   on   a    specific   atom 

FA:  The   component    force    increment    on   an   atom 

FDTI:    DTI    X  CVD,    a   parameter    used   to  determine   DT   by   maximum 

energy   method 
FM:  A    small    number    used    in   checking  potential    energy    zero 

point 
FM2:  FM    squared 

FMAX:  Maximum   total    force    on    the  most    stressed  atom    in    the 

crystal 
FMAX1 :         Maximum   total    force   on   Atom   1 

FOD:  FORCE/DIST--a    ratio   used   to   avoid  repeated   division 

FORCE:         Numerical    value   of   the    force   function   with  a    variable 

parameter 
FX(I),    FY(I),    FZ(I):    x,y,z    components    of   total    foce    on   an   atom 
FXA :  Born-Mayer    force   function  parameter 

HBMAS:  ^  BMAS--  a  ratio  used  to  avoid  repeated  division 
HDTOD:  \  DTOD--  a  ratio  used  to  avoid  repeated  division 
HDTOM:  \  DTOM--  a  ratio  used  to  avoid  repeated  division 
HDTOMB:  ^  DTOMB--  a  ratio  used  to  avoid  repeated  division 
HTMAS:  \  TMAS--  a  ratio  used  to  avoid  repeated  division 
II:  Variable    in    cubic    fit    subroutine 

13:  Variable    in    cubic    fit    subroutine 

I  DEEP:         Number    of   mobile    layers 
IH1 :  Alpha    numeric   array    for   program   title 

IH2 :  Alpha    numeric   array    for    Morse   function   parameters 

IHB:  Alpha    numeric   array    for    bullet    element 

IHS:  Alpha   numeric   array    for    type  and   orientation    of   crystal 
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IHT:      Alpha  numeric  array  for  target  element 

ILAY:     Same  as  I DEEP 

IN:       Odd-even  integer  used  to  determine  atom  site  establishment 

IP:       Subscript  value  of  atom.   Used  in  subroutine  STEP  and 

ENERGY 
IQ:       Parameter  that  determines  whether  or  not  a  self  defect  is 

to  be  given  a  repulsive  potential  or  a  composite  attractive- 
repulsive  potential 
ISHUT:    A  parameter  used  to  shut  down  the  program 

IT:       Unsealed  fixed  point  x  coordinate  used  in  lattice  generation 
ITT:      Odd-even  integer  used  to  determine  atom  site  establishment 
ITYPE:    Parameter  used  to  determine  the  type  of  point  defect: 

vacancy,  self-interstitial,  replacement,  foreign  inter- 
stitial 
IVACX,  IVACY,  IVACZ:  Input  plane  numbers  to  specify  NVAC 
IX,  IY,  IZ:  Number  of  x,y,z  planes  of  crystal 
J2 :       Variable  in  the  cubic  fit  subroutine 
JT:       Unsealed  y  coordinate  used  in  crystal  generation 
JTS:      Variable  used  to  establish  atom  sites 
JTT:      Variable  used  to  establish  atom  sites 
KF:       Final  K  in  LOCAT  (K)  assigned  to  an  atom 
KT:       Unsealed  z  coordinate  used  to  establish  atom  site 
LCUT(I):  Used  to  identify  an  ith  atom  which  is  not  included  in 

calculations 
LD:       The  highest  numbered  atom  in  the  mobile  layers 
LL:       The  highest  numbered  atom  in  the  entire  crystal 
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LOCAT(K) : Dimensioned  variable  that  remembers  the  numbers  of  the 
atoms  within  a  radius  ROEL  of  the  primary  at  time  zero 

LS:  Variable  associated  with  each  of  the  nine  lattice  gene- 
rator subroutines 

MCRO:  One  number  higher  than  the  order  of  the  fit  between  the 
Born-Mayer  and  Morse  potentials,  always  4  in  this  simu- 
lation 

ND:       Data  output  increment,  in  numbers  of  timesteps 

NDEC:     Counting  index  for  DTI  variation 

NEW:  Parameter  used  to  determine  whether  or  not  atom  numbers 
have  been  stored  in  LOCAT(K) 

NPAGE:    Page  numbering  variable 

NRUN:     Parameter  used  to  determine  whether  or  not  to  read 
additional  data  cards 

NS:       Initial  print  statement  timestep  number 

NT:       Timestep  number 

NTT:      Timestep  number  limit  before  shutdown 

NVAC:     An  atom  number  used  to  establish  point  defects  or  used  as 
a  reference  point  for  interstitial  placement 

PAC:      Parameter  for  bullet  force  function  correction 

PBMAS:    Primary  mass  in  kilograms 

PEXA,  PEXB:  Input  Born-Mayer  potential  function  parameters  for  the 
bullet -target  interaction 

PFPTC:    Primary  force  function  evaluated  at  ROE 

PFXA:     Primary  force  function  parameter 

PKE(I):   Kinetic  energy  of  the  ith^  atom 

PLANE:    Alpha -numer ic  array  for  lattice  orientation 
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POT:      Potential  energy  between  two  atoms 

PPE(I):   Potential  energy  of  the  ith  atom 

PPTC:     Primary  potential  function  evaluated  at  ROE 

PTE(I):   Total  energy  of  the  ith  atom  (potential  +  kinetic) 

PTMAS:    Target  mass  in  kilograms 

RE:       Input  Morse  potential  parameter 

RLL:      Reciprocal  of  LL 

RO:       Spacing  constant  in  FCC(llO)  lattice  generation  subroutine 

ROE:      Nearest  neighbor  distance 

R0E2:     ROE  squared 

ROEA:     Maximum  cut  off  for  Born -Mayer  potential 

ROEB:     Minimum  cut  off  for  Morse  potential 

ROEC:     Maximum  cut  off  for  Morse  potential 

ROEC2:    ROEC  squared 

ROEL:     Radius  inside  of  which  local  potential  energy  is  found 

R0EL2:    ROEL  squared 

ROEM:     ROE-DTI,  region  in  which  modification  of  repulsive  force 

must  be  made 
RX(I),  RY(I),  RZ(I):  x,y,z  coordinates  of  an  ith  atom  at  any  time 
RXI(I),  RYI(I),  RZI(I):  x,y,z  coordinates  of  an  ith_  atom's  initial 

position 
RXK(I),  RYK(I),  RZK(I):  x,y,z  coordinates  of  temporary  position  of 

an  ith_  atom  during  force  cycle 
SAVE:     h   POT 

SCX,  SCY,  SCZ:  x,y,z  coordinate  scale  factors 

START:    An  optional  timing  variable,  not  used  in  this  simulation 
SUM:      Variable  in  cubic  fit  subroutine 
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TARGET:   Alpha -numer ic  array  for  target  material 

TSAVE:    Bullet  mass/(target  mass  +  bullet  mass);  distributes 

potential  energy  between  target  and  bullet 
TE:       Total  energy  of  all  crystal  atoms  (kinetic  +  potential) 
TEMP:     Temperature  of  lattice  in  degrees  Kelvin.   Not  used  in 

this  simulation 
TEFAC:    Product  of  DTI  *  CVD 
TFAC:     A  time  factor  ratio  used  to  determine  DT  by  maximum  force 

method 
TFACB:    TFAC  for  the  bullet 

THERM:    Thermal  energy  of  atom.   Not  used  in  this  simulation 
TIME:     Elapsed  problem  time  in  seconds 

TLPE:     Total  local  potential  energy  of  atoms  within  a  radius  ROEL 
TLPE0:    TLPE  at  time  zero 
TMAS:     Target  atom  mass  in  amu 

TPKE:     Total  kinetic  energy  of  all  crystal  atoms 
TPOT:     Total  potential  energy  of  all  crystal  atoms 
TPCTL:    Storage  position  for  the  last  computed  value  TPOT 
VSS :      Storage  variable  for  velocity  components 
VX(I),  VY(I),  VZ(I):  x,y,z  components  of  ith_  atoms  velocity 
X,  Y,  Z:  Unsealed  coordinates  used  in  crystal  generation 
XSTART:   X  Coordinate  used  in  impact  point  generator 
YLAX(I):  Relaxation  in  -y  direction  of  ith_  layer  in  L.U. 
ZP:       Floating  point  form  of  JTT 
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COMPUTER  PROGRAM 

THE  FOLLOWING  PROGRAM  IS  THE  ONE  USED  FOR  THE  STATIC 
SIMULATION  RUNS.   THE  DYNAMIC  PROGRAM  IS  ESSENTIALLY  THE 
SAME  AS  THE  STATIC.   IN  THE  DYNAMIC  PROGRAM,  SUBROUTINES 
(BlOO)  AND  (PLACE)  ARE  OMITTED,  AND  THE  INITIAL  ATOM  POSI- 
TIONS ARE  READ  IN  ON  DATA  CARDS.   ALSO, AN  IMPACT  POINT 
GENERATOR  PACKAGE  IS  INCLUDED  IN  THE  DYNAMIC  PROGRAM. 
BOTH  PROGRAMS  CALCULATE  THE  EQUATIONS  OF  MOTION  FOR  THE 
SYSTEM  THROUGH  THE  USE  OF  THE  APPROPRIATE  POTENTIAL  FUNCTION 
AND  THE  AVERAGE  FORCE  METHOD  OF  INTEGRATION.   THE  NET 
DISTANCE  OF  MOVEMENT,  VELOCITY,  AND  ENERGY  VALUES  ARE 
PRINTED  OUT  FOR  THE  DESIRED  ATOMS  ON  SELECTED  TIMESTEP 
NUMBERS.   THE  FINAL  POSITION,  VELOCITY,  AND  ENERGY  VALUES 
FOR  EACH  ATOM  ARE  SUMMARIZED  AT  THE  END  OF  THE  PROGRAM. 
//  EXEC  FORTHCLG, TIME. GO=20, REGION. GC=140K 
//FORT.SYSIN  DD  * 
C 

DIMENSION  VX(IOOO) , VY( 1000 ) ,VZ ( 1 OOO ), PKE( 1 OOO ) 
DIMENSIONING  OF  VARIABLES  NOT  NEEDED  IN  COMMON 
C 

DIMENSION  DX(IOOO) ,DY( 1000) , DZ ( 1 000 ) , P TE ( 1 000 ) 

DIMENSION  RXK(IOOO) , RYK ( 1000 ) , RZK ( 1000} 
COMMON  LABELING  OF  VARIABLES  REQUIRED  IN  OTHER  SUBROUTINES 

C0MM0N/C0M1/RX( 100  0),RY(1000),RZ(1000),LCUT(  1000)  , 
1LL,LD,ITYPE,NVAC 

COMMON /COM  2/IHK  2  0  )  ,  I  H2(  8  )  ,  I  HS  (  1  0)  ,  I  HB  (  6)  ,  I  HT  {  6  )  , 
1TARGET (4) , TM AS , BULL ETC  4) , BH AS , PL AN E , TEMP  ,  THERM 

COMMON /C0M3/RX I ( 1 0  00 ) , RY I { 1 000 ) , RZ I ( 10  00 ) , CVR, EVR, 
INT, TIME, DT, DTI , ILAY 
1IVACX, IVACY, IVACZ 

COMMON /C0M4/ IX,IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z, 

COMMON /C0M5/R0E,R0E2,R0EM,EXA,EXB, PEXA,PEX8,FXA,PFXA, 
1 IQ,TSAVE,BSAVE 

COMMON /C0M6/FXC 100  0) ,FY(1000) , FZ(IOOO) ,PAC,PFPTC, FM 

COMMON/ C0M7/P PTC, T POT, PPE( 1 000 ) , TLP E , R OEL , RC5L2 ,NE W 

COMMON/COM8/ROEA,ROEB,ROEC,ROEC2,CPO,CP1,CP2,CP3, 
1CF0,CF1,CF2,CGD1,CGD2,CGB1 ,CGB2 ,CGF1 ,CGF2 

COMMON/COMA/  A(4,5),MCP0 
C 
READ  STATEMENT  FORMATS 

9010  FORMAT ( 20 A4) 

9020  FQRMAT( 3A4,3FS. 5,2^5.2) 

9030  F0RMAT(4A4,3F8.5,6A4,F6.2) 

9040  FORMAT (F6. 3, 5X, 15, 614, 3 F5, 312) 

90  50  FORMAT! 10A4, A4 , 41 3 , F 8. 4, 14 , F5. 3 ) 
C 

WRITE  STATEMENT  FORMATS 

9610  FORMAT (1H1) 

9620  FQRMAT(47X,«  SUMMARY  OF  ATOMS '//, 35 X , 8A4, « ,  NT  =  M4,//, 

13(  '  ATOM       POSITION       BIND  ENERGY    '),//) 
9  630  FORMAT (3(I5,3F6.2,F3.4,8X)) 
9640  FORMAT (/4X,F10. 3, 25H  EV, TOTAL  KINETIC  ENERGY ,, F10 . 3, 

1H  EV, TOTAL  POTENTIAL   ENERGY ,F 1 0. 3 , '  EV,  REDUCTION', 

1//60X'  RADIUS  =  '  ,r-5.2,  ) 
96  50  FORMAT ( 105X,4HPAGE,I3,/,1H1) 
9660  FORMAT (/   •      ATOM     DX         DY         DZ 

1VX         VY         VZ        KE         PE         TE'/) 
9670  FORMAT! 1 1  8 ,3F 1 0. 3 , 3F 10. 1,3F10.4  ) 
9680  FORMAT (•  SHARP  DT  DECREASE ', 2E 10.3 ) 

9690  F0RMAT(I4,3F5.2,I4) 

9691  FORMAT! 9F8. 4) 

9692  FORMAT  (IX,  14,/) 
C 

INITIALIZING 

START=0.01*ITIME(XX) 

DO  2  1=1,1000 
RXK(I)=0.0 
RYK( I ) =0.0 
RZK( I ) =0.0 
VX( I )=0.0 
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c 

INPUT  DATA 


VY(I )=0.0 
VZ( I  1  =  0.0 

PKE( 11=0.0 
PPE(  11=0.0 
PTE{ II =0.0 
RZK  11=0.0 
ISHUT=1 
NRUN=0 


READ  (  5,9010)  IH1 
READ  (  5,90201  IH2 


L 
.2,DC0N,ALPHA,RE,R0EC,R0EL 
READ  (  5,9030)  BU LLET , BMAS , PEXA ,P EXB, I H8, THERM 

*1GET,TMAS,EXA,EXB,IHT,TEMP 

;.PI  AMP. I  q  .  T  Y  .  T  Y  .  T  7  .r\/R  .  MTI 


READ  (  5,9030)  TARGET , TMA S ,EXA , EXB , I HT , TEMP 

READ  (  5,9050)  IHS, PL ANE, L S , I  X, I Y , I Z , CVR , MCRC  ,DTI 

C 

CONSTANTS  AND  SCALING  FACTORS 

DTIS=DTI 

RCE2=4.0 

R0E=SQRT(R0E2) 

ROEM  =  ROE-DTI 

R0EL2=R0EL*R0EL 

CVE=1.60E-19 

CVM=1.672E-27 

VFAC=.50 

FM=1.0E-10 

FM2=FM*FM 

CVD=CVR*1.0E-10 

CVED=CVE/CVD 

PTMAS=TMAS*CVM 

PBMAS=BMAS*CVM 

HTMAS=0.5*PTMAS/CVE 

HBMAS=0.5*PBMAS/CVE 

TSAVE=BMA5/( BMAS+TMAS) 

BSAVE=TMAS/(BMAS+TMASJ 
C 
REPULSIVE    POTENTIAL    PARAMETERS 

FXA=ALOG(-EXB*CVED)+EXA 
PFXA=A LOG {-PEXB*CVED)+ PEXA 

PPT C= EXP (PEXA+PEXB*ROE) 
PAC=ALOG(CVED)+PEXA 

PFPTC=EXP(PAC+PEX3*R0E) 
C 
ATTRACTIVE    POTENTIAL    PARAMETERS 

CGD1=AL0G(DC0N1+2.0*ALPHA*RE 
CGD2=AL0G(2.0*DC0N)+ALPHA*RE 
CGB1=-2.0*ALPHA*CVR 

CGB2=-ALPHA*CVR 
CGF1=AL0G<-CG31*CVED1+CGD1 
CGF2=AL0G(-CGB2*CVED1+CGD2 
C 

CUTOFF    DISTANCES    FOR    ATTRACTIVE    AND    REPULSIVE    POTENTIALS 

R0EA=1.50/CVR 

R0EB=2.0/CVR 

R0EC2=R0EC*R0EC 
C 

PARAMETERS    FOR    CALCULATION    OF    THE    BEST    CUBIC    FIT     IN    THE    GAP 
BETWEEN    MAXIMUM    DISTANCE    CUTOFF    OF    THE    REPULSIVE    POTENTIAL 
(ROEA),     AMD    MINIMUM    DISTANCE    CUTOFF    OF    THE    ATTRACTIVE    POTEN- 
TIAL    (ROEBJ.       SUBROUTINE    CROSYM    ACTUALLY    PERFORMS    THIS    CURVE 
FITTING. 

A(l, 11=1.0 

A(1,21=R0EA 

A(  1,3)=R0EA*R0EA 

A(1,4)=RGEA"*3 

A(1,51=EXP(EXA+EXB-RCEA) 

A(2, 11=1.0 
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A(2,2)=ROEB 
A(2,3)=R0EB*R0EB 
A(2»4)=R0EB*#3 

A(2  ,5)=EXP(CGD1+CG31*R0EB)-EXP(CGD2+CGB2*R0EB) 
A(3,l)=0.0 
A(3,2)=-1.0 
A(3,3)=-2.0*R0EA 
A( 3,4) =-3. 0*ROEA*ROEA 
A(3, 5)=EXP(FXA+EXB*R0EA)/CVED 
A(4,l)=0.0 
A(4,2) =-1.0 
A(4,3)=-2.0*R0E3 
A(4,4) =-3.0*R0EB*R0EB 

A(4,5) =(EXP(CGF1+CGB1*R0EB)-EXP(CGF2+CGB2*R0EB) )/CVED 
CALL  CROSYM 
CPO=A( 1,5) 
CP1=A(2,5) 
CP2=A(3,5) 
CP3=A(4,5) 
CF0=-CP1*CVED 
CF1=-2.0*CP2*CVED 
CF2=-3.0*CP3*CVED 
5   READ  (  5,9040)  EVR , NTT,NS ,ND , I P , I  DEEP , I  TYPE , NVAC , D1X, 
1D1Y,D1Z, IVACX, IVACY, IVACZ 
IF(NTT.EQ.O)  GO  TO  9999 
IQ=ITYPE-1 
EV=EVR*1.0E+3 

SELECTION  OF  THE  DESIRED  CRYSTAL  STRUCTURE  AND  ORIENTATION. 
100,  110,  AND  111  PLANES  OF  FACE-CENTERED,  BODY-CENTERED, 
AND  DIAMOND  STRUCTURES  ARE  ALLOWED.   ILAY  AND  IDEEP  ARE  VAR- 
IABLES ESTABLISHING  THE  NUMBER  OF  MOBILE  LAYERS  IN  THE 
CRYSTAL.   RXI(I)  AND  RXK  (  I  )  ARE  VARIABLES  SAVING  THE  ORIGIN- 
AL X-POSITION  OF  THE  I ■ TH  ATOM.   Y  AND  Z  POSITIONS  ARE 
ANALOGOUS. 

14  CALL  B100 
30  ILAY=IDEEP 

IF(IDEEP)  35,35,40 
35  LD  =  LL 

ILAY=IY 
40  RLL=1.0/LL 
TP0TL=1.0 
DO  45  1=1, LL 
RXK( I)=RX( I) 
RYK( I ) =RY(I) 
RZK(I)=RZ(  I) 
RXI ( I J=RX(I) 
RYI ( I ) =  RY( I) 
45  RZI ( I )=RZ(  I) 
C 

THIS  SECTION  ALLOWS  ONE  TO  REPEAT  A  RUN  OF  THE  PROGRAM  WITH 
DIFFERENT  DATA  WITHOUT  REPEATING  INITIALIZATION,  POTENTIAL 
PARAMETER  CALCULATIONS  AND  CRYSTAL  LATTICE  BUILDING.   SUB- 
ROUTINE PLACE  USES  LCUT(I)  AND  NVAC  TO  CREATE  VACANCIES, 
INTERSTITIALS,  AND  REPLACEMENT  IMPURITIES  AT  DESIRED  LOCA- 
TIONS IN  THE  LATTICE. 

IF(NRUN.EO.O)  GO  TO  60 

DO  55  1=1, LL 

LCUTU  )=0 

RX(I)=RXI( I) 

RY( I)=RYI ( I) 

RZ( I ) =  R Z I (  I) 

RXK(  I )=RXI (I  ) 

RYK(  I)=RYI  II) 
55  RZK( I)=RZI (I) 
60  NRUN=1 

CALL  PLACE 

RXK  1)  =RX(  1) 

RYI (1)=RY<  1) 

RZI (1) =RZ(1) 
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RXK( 1)=RX( 1) 

RYK( 1)=RY(1) 

RZK(l) =RZ(1) 

DO  65  1=1, LL 

VX(I )=0.0 

VY(I }=0.0 

VZ(I )=0.0 

PPE( I  J =0.0 

PKE( I) =0.0 
65  PTE( I)=0.0 

TPOT=0.0 

NEW  =  0 
C 

THE  ENERGY  SUBROUTINE  CALCULATES  THE  POTENTIAL  ENERGY  OF 
EACH  ATOM  IN  THE  LATTICE.   SUBROUTINE  LOCAL  SUMS  UP  THIS 
ENERGY  FOR  ALL  ATOMS  WITHIN  A  SPECIFIC  RADIUS  OF  THE  POINT 
DEFECT. 

THIS  SECTION  PRINTS  OUT  X,  Y,  AND  Z  COORDINATES,  IN  LATTICE 
UNITS,  AND  BINDING  ENERGIES  OF  EACH  ATOM  IN  THE  CRYSTAL  AT 
TIME  ZERO. 

CALL  ENERGY 

BIND=-TPOT 

TE=TPOT+BIND 
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C 

THIS  I 

FORCE 

THE  DY 

CALCUL 

FORCES 

THE  PR 

STEP; 

AND  AL 

FORCES 

CULATI 

CALCUL 

STEP, 

ALLOWE 

TIMEST 


TIME  = 
NT  =  0 

WRIT 

WRIT 
DO  70 
K=I  +  1 
J=I  +  2 

WRIT 
RY(K) 

WRIT 
NPAGE 
NPAGE 

WRIT 


0.0 

E  (  6,9610) 

E  (  6,9620) 

1=1, LL, 3 


IH2,NT 


E  (  6,9630)  I ,RX( I)  ,RY( I) ,RZ(  I)  ,PPE( I  )  ,K,RX(K), 
i  RZ(KJ*PPE(K)  *JrRX(  J)  rRYt  J)  tRZ(  J)  tPPEU* 

E  (  6,9640)  TPKE,TPOT,TE,ROEL 

=  1 

=NPAGE+1 

E  (  6,9650)  NPAGE 
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OF  EACH 


DT=1.0E-15 

DDTI=.005 

NDEC=0 
100    DTI=DTI-DDTI 

NDEC=NDEC+1 

DTOD=DT/CVD 

TFAC=2.0*PTMAS*DTI*CV0 

TFACB=2.0*PBMAS*DTI*CVD 

TEFAC=DTI*CVD 

HOTOD=0.5*DTOD 

DTOM=DT/PTMAS 

HDTOM=0.5*DTOM 

DT0M3=DT/PBMAS 

HDT0MB=0.5*DTCMB 
200  CALL  STEP 

IFUCUT(l)  .GT.O)  GO  TO  240 

1  =  1 

RXK( I }=RX( I) 
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240 


245 


260 
265 


275 


280 


300 
310 
320 


RYK( I) 
RZK(I) 
RX(I )  = 
RY  (  I  )  = 
RZ(I )= 
DO  245 
IFUCU 
RXK( I ) 
RYK( I ) 
RZK( I) 
RX(I )= 
RY ( I  J  = 
RZ(I  )  = 
CONTIN 
CALL  S 
EMAX=0 
FMAX=0 
TIME=T 
NT=NT+ 
IF(LCU 
1=1 
VSS=VX 

vx(  n  = 

RX(  I  )  = 
VSS=VY 
VY( I )= 
RY(  I  )  = 

vss=vz 

VZ  (  I  )  = 
RZ(  I  )  = 
PKE( I ) 
EMAX1= 
FMAX1= 
DTE1=T 
DTF1=S 
FX(  I  )  = 
FY( I )  = 
FZ(  I)  = 
DO  280 
IF(LCU 
VSS=VX 
VX(I )  = 
RX( 1 ) = 
VSS=VY 
VY(I )  = 
RY( I )= 

vss=vz 

VZf 1  )  = 

RZ( I)= 
PKF( I ) 
F2=FX( 
FX(  I  )  = 
FY (I  )  = 
FZ(  I  )  = 
IF(F2. 
IF(PKE 
CONTIN 
DTL=DT 
CTIME= 
DTE=TE 
DTF=SQ 
IF(EMA 
DT=DTE 
IF(DT. 
IF(DT. 
IF(DT. 
IF(DT. 
IF(  ISH 
IF  (IMS- 
DO  325 
VX(  I)  = 
VY(I)= 


=  RY 
=  RZ 
RX( 
RY( 
RZ( 
1  = 
T(I 
=  RX 
=  RY 
=  RZ 
RX( 
RY( 
RZ( 
UE 
TEP 
.0 
.0 
I  ME 
1 
T(l 


(I) 
(I) 
I  )  + 
I)  + 
I  )  + 
2,L 
).G 
(IJ 
(  I) 
(  I) 
I)  + 
I  )  + 
I)  + 


DTCD*{HDT0M8*FX( I)+VX( I ) ) 

DTGD*(HDTOMB*FY<  I  )+VY(  I )  ) 

DTQD*(HDTOMB*FZ( I)+VZ( I ) ) 

D 

T.OGO    TO    245 


DTOD*(HDTOM*FX( I )+VX( I ) ) 
DTOD*( HDTOM*FY( I  )  +  VY(I  ) ) 
DTOD*(HDTOM*FZ(  I)+VZ( I ) ) 


+  DT 

)  .GT.OJ     GO    TO    265 


(I) 

VSS+HDTO 
RXK( I)+( 
(I  ) 

VSS+HDTO 
RYK( I )+( 
(I) 

VSS+HDTO 
RZK( I)+{ 
=VX( I)*V 
PKE( I) 
FX(I)*FX 
EFAC-SGR 
QRT(TFAC 
0.0 
0.0 
0.0 

1=2 
TCI) 
(I  ) 

VSS+HDTO 
RXK(  I)4( 
(  I  ) 

VSS+HDTO 
RYK(I)+( 
(I) 

VSS+HPTO 
RZK(I)+( 
=VX( I)*V 
I)*FX( I ) 
0.0 
0.0 
0.0 

GT.FMAX) 
(  I  )  .  GT  .  E 
UE 


iLD 

.GT. 


MB*FX( I ) 

VX( I J+VSS)*HDTOD 

MB*FY( I ) 

VY(I )+VSS)*HDTCD 

MB*FZ< I ) 

VZ( I  J+VSS  )*HDTOD 

X(I)+VY( IJ*VY(I )+VZ(  I  )*VZ( I) 

(I )  +  FY( I )*FY(  I  )+FZ(I  )*FZ( I ) 

T( 1.0/FMAX1) 

B/FMAX1J 


0)G0    TO    280 

M*FX( I ) 

VX(I )+VSS)*HDTCD 

M*FY( I ) 

VY(I  J  +  VSSKHDTOD 

M*FZ(i) 

VZ( I )+VSS)*HDTOD 

X( I )+VY( I )*VY(I )+VZ( I )*VZ(I ) 

+  FY(  I )*FY( I)+FZ( I  )*FZ(  I) 


FMAX=F2 

MAX)     EMAX=PKE(I) 


0.01*ITIME(XX)-START 
FAC*SQRT(1 .O/EMAX) 
RTtTFAC/FMAX) 
Xl.GT.EMAX)     EMAX=EMAX1 
1 

GT.DTF1)     DT=DTF1 
GT.DTE)     DT=DTE 
GT.DTF)     DT=DTF 
GT.1.0E-14)     DT=1.0E-14 
UT.EQ.-l)     GO    TO    400 
NT)     400,400,320 

1=1, LL 
VFAOVX(  I) 
VFAC*VY( I ) 
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325  VZ( I )=VFAC*VZ( I ) 
GO  TO  800 
C 

THE  PRINT  SUBROUTINE  PLACES  A  HEADING  OF  PERTINENT  INFORMA- 
TION AT  THE  TOP  OF  EACH  TIMESTEP  PRINTOUT. 

POTENTIAL  ENERGY  AND  LOCAL  POTENTIAL  ENERGY  FOR  EACH  ATOM 
ARE  CALCULATED  BASED  ON  THE  NEW  POSITIONS.   SUMMATIONS  OF 
TOTAL  POTENTIAL  AND  KINETIC  ENERGY  FOR  THE  LATTICE  ARE  PER- 
FORMED.  DX,  DYt  AND  DZ  KEEP  TRACK  OF  MOTION  RELATIVE  TO  THE 
INITIAL  POSITION  AT  TIME  ZERO  FOR  EACH  ATOM. 

400  CALL  PRINT 


410 
450 


620 


700 


TPOT=0 
DO  450 
PPE( I ) 
PTE(  I) 
CALL  E 
PKE(l) 
TPKE=P 
PTE( 1) 
DO  620 
PKE( I) 
TPKE=T 
PTE(  I  ) 
TE=TPO 
WRITE 
DTEST= 
IF  (DT 
IF1TP0 
ERAT=T 
DO  750 
OX (I  )  = 
DY  (  I  )  = 
DZ(  I  )  = 
IF  (DX 
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ENERG 
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IF 
IF 
GO 
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TO 


.0 

1=1 
=  0.0 
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NERG 
=HBM 
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1=2 
=  HTM 
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=  PKE 
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RX(  I 
RY(  I 
RZ(  I 
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ND 

t9 

1) 
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*PKE( 1) 

)+PPE( 1) 
L 

*PKE(I  ) 
Ed  ) 
)+PPE( I) 


660) 
-RYI 
.  0. 
POTL 
L 
D 

RXI  ( 
RYI  ( 
RZI  ( 
.GE. 
.GE. 
.GE. 


(1)  )**2 

01)  DTEST=  0.01 
)  GO  TO  700 


I) 

I) 
I  ) 

DTEST)  GO  TO  720 
DTEST)  GO  TO  720 
DTEST)  GO  TO  720 


SECTION  PRINTS  THE  RELATIVE  MOTION,  VELOCITY,  AND 
Y  OF  EACH  ATOM,  FOR  EVERY  TIMESTEP  SO  DESIGNATED:   IE 
ND'TH  TIMESTEP,  BEGINNING  WITH  #NS  AND  ENDING  WITH 


720 
750 


760 


780 

790 
800 

810 


950 
C 

THIS  S 
ENERGI 
PROGRA 
OF  EAC 


WRITE 
VZ(  I  )  , 
CO NT  IN 
WRITE 
WRITE 
NPAGE= 
TPOTL= 
IF(NT- 
DO  780 
VX(I )= 
VY(I  )  = 
VZ(I  )  = 
CONTIN 
IF( ISH 
NS=NS+ 
IF  (ND 
GO  TO 
DDTI=0 
DTI=DT 
NDEC=0 
GO  TO 
CONTIN 


( 

PK 

UE 

( 

( 

NP 
TP 
NT 
I 
VF 
VF 
VF 
UE 
UT 
ND 
tC 
10 
.1 
1  + 


6,967  0)     I ,DX( I) ,DY( I),DZ( I ),VX( I) ,VY( I), 
E( I) ,PPE(I ) ,PTE( I) 

6,9640)     TPKE,TPOT,TE,ROEL 

6,9650)     NPAGE 
AGE  +  1 
OT 

T)  760,950,950 
=  1,LL 
AC-VXU  ) 
AC*VY(I ) 
AC*VZ( I ) 

.EQ.-l)  GO  TO  950 

.EQ.10)  GO  10  810 
0 

*DDTI 
DDTI 


100 
UE 


ECTION  PRINTS  OUT  X,  Y,  AND  Z  COORDINATES  AND  BINDING 
ES  OF  EACH  ATOM  IN  THE  CRYSTAL  AT  THE  END  OF  THE 
M.   ALSO, DATA  CAPOS  ARE  PRINTED  WITH  X,Y,Z  COORDINATE 
H  ATOM  IN  THE  CRYSTAL  FOR  USE  IN  THE  DYNAMIC  PROGRAM. 
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955   WRITE  (  6,9620)  IH2,NT 

WRITE  (7t9690)  LL , D1X , D1Y , D1Z , NVAC 
DO  965  1=1, LL, 3 
K=I  +  1 
J=I+2 
WRITE  (7,9691)  PX( I ) , RY( I ) , RZ ( I ) , RX ( K) , RY ( K ) , RZ ( K ) , RX 
1RZ(J ) 
965   WRITE  (  6,9630)  I , RX (  I  ) , RY ( I  )  , RZ(  I  )  , PPE ( I  )  , K ,RX( K ) , 
1RY(K) ,RZ(K),PPE(K) ,J,RX(J) ,RY(J),RZ(J),PPE(J) 
WRITE  (  6,9640)  TPKE , TPOT , TE, ROEL 
WRITE  (  6,9650)  NPAGE 
1000  IF(ISHUT)  9999,5,5 
9999  STOP 
END 

SUBROUTINE  CROSYM 
C 

SOLVES  M  SIMULTANEOUS  EQUATIONS  BY  THE  METHOD  OF  CROUT 
THIS  SUBROUTINE  FITS  THE  BEST  CUBIC  BETWEEN  THE  REPULSIVE 
AND  ATTRACTIVE  PARTS  OF  THE  POTENTIAL. 
C 

COMMON/COMA/  A(4,5),MCR0 

M=MCRO 

N=M+1 

11  =  1 
100  13=11 

SUM=ABS(A( 11,11)) 

DO  120  1=11, M 

IF(SUM-ABS(A(I ,11)))  110,120,120 
110  13=1 

SUM=ABS( A( 1,11)) 
120  CONTINUE 

IF(13-U)  130,150,130 
130  DO  140  J=1,N 

SUM=-A( II, J) 

A(  II, J  )  =  A(  13, J) 
140    A( I3,J)=SUM 
150    13=11+1 

DO  160  1=13, M 
160  A( 1,11 )=A( I, II )/A( 11,11) 
170  J2=I 1-1 

13=11+1 

IF(J2)  180,200,180 
180  DO  190  J=I3,N 

DO  190  1  =  1,  J2 
190  A(  II , J)=A(  II ,J)-A{  II, I  )*A( I , J) 

IF(Il-M)  200,220,200 
200  J2=I1 

11=11+1 

DO  210  1=1 1,M 

DO  210  J=1,J2 
210  A(  I  ,I1)=A( I, Il)-Al I ,J)*A(J,  II) 

IF(Il-M)  100,170,100 
220  DO  240  1=1, M 

J2=M-I 

I3=J2+1 

A(  I3,N)  =  A( I3,N)/A(  13, 13) 

IFU2)     230,250,230 
230    DO    240    J=1,J2 

240    A( J,N)=At J,N)-A(  I3,N)*A{ J,  13) 
250    RETURN 

END 
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SUBROUTINE  BlOO 
THIS  IS  A  LATTICE  GENERATOR  THE  THE  BCC  (100)  ORIENTATION. 
THE  CRYSTAL  IS  DEVELOPED  IN  THE  ORDER,  X  FOLLOWED  BY  Z, 
FOLLOWED  BY  Y. 

IT  CONTAINS  A  NONSTANDARD  USE  OF  THE  SURFACE  RELAXATION 
PARAMETER. 
C 

COMMON/CCM1/RX(1000),RY( 1000 ), RZ ( 1000 ), LCUT ( 1000) , 
1LL,LD, ITYPE,NVAC 

C0MM0N/C0M4/IX, IY,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z, 
HVACXt  IVACY,  IVACZ 

DIMENSION  YLAX(20) 

DATA    YLAX/20*0.0/ 

YLAXd  )=-0.20 

YLAX(2)=-0.03 

SCX-1.0 

SCY=1. 0 

SCZ=1.0 

M    =    2 

JT=0 

Y=-SCY 

DO    60    J=l, IY 

Y=Y+SCY 

KT  =  0 

Z=-SCZ 

DO    59    K=1,IZ 

Z=Z+SCZ 

IT=0 

X=-SCX 

DO    5  8    1  =  1,  IX 

x=x+scx 

IF( IT-( IT/2)*2)     21,11,21 

11  IF(JT-(JT/2)*2)     57,12,57 

12  IF(KT-(KT/2)*2)     57,30,57 

21  IF(JT-( JT/2)*2)     22,57,22 

22  IF(KT-(KT/2)-2)     30,57,30 
30    RX(M)=X 

RY(M)=Y+YLAX(J) 

RZ(M)=Z 

M  =  M+1 

IF    ( IT.NE.IVACX)     GO    TO    57 

IF     ( JT  .NE.  IVACY)    GO    TO    57 

IF     (KT.NE. IVACZ)     GO    TO    57 

NVAC=M-1 

57  IT=IT+1 

58  CONTINUE 
KT=KT+1 

59  CONTINUE 

JT    =    JT    +    1 
IF(IDEEP-JT)     60,110,60 

60  CONTINUE 
LL=M-1 

100    RETURN 
110    LD=M-1 

GO    TO    60 

END 

SUBROUTINE  PLACE 
C 

THIS  SUBROUTINE  LOCATES  A  VACANCY,  INTERSTITIAL,  OR  REPLACE- 
MENT IMPURITY  IN  THE  LATTICE. 
C 

COMMON/COM 1/RX( 100  0),RY(100  0),RZ(1000),LCUT(10  00), 
1LL,LD,  ITYPE,MVAC 

C0MM0N/CCM4/IX,I  Y ,  I  Z  ,  SCX  ,  SCY  ,-SCZ  ,  I  DEEP  ,  Dl  X  ,  D1Y  ,  Dl  Z  , 
1IVACX, IVACY, IVACZ 

GO  TO   (10,20,30,40),  ITYPE 
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10  LCUT(NVAC)  =  1 

LCUTU)  =  1 

RX(1)=0.0 

RY(1)=0.0 

RZ(1)=0.0 

GG  TO  50 
20  RX(1)=RX(NVAC)+D1X 

RY(1 )=RY(NVAC)+D1Y 

RZ(1)  =  RZ(NVAC)  +  D1Z 

GO  TO  50 
30  LGUT(NVAC)  =  1 

PX(1)  =  PX(NVAC) 

RY(1)  =  P.YINVAC) 

RZ(1)  =  RZ(NVAC) 

GO  TO  50 
40  RX(1)=RX(NVAC)+D1X 

RY(1)=RY(NVAC)+D1Y 

RZ(1)  =  RZ(NVAC)+D1Z 
50  RETURN 

END 

SUBROUTINE  STEP 
C 

THIS  SUBROUTINE  DOES  THE  DYNAMICS  FOR  ONE  TIMESTEP. 
THE  FIRST  HALF  DOES  THE  DYNAMICS  FOR  ATOM  #1;  THE  SECOND 
HALF  FOR  ALL  OTHERS. 
C 

COMMON/COM 1/RX( 1000) ,RY( 1000 ), RZ ( 1 000 ) ,LCUT( 1000) , 
1LL,LD, ITYPE,NVAC 

COMM ON /CCM5/R0E, RO E2, ROEM , E XA, E XB,P EX A, PEXB,FX A, PFX A, 
1IQ,TSAVE»BSAVE 
COMMON/CO,'i6/FX(100  0)  , FY(1000),FZ(1000),PAC,PFPTC,FM 
C0MMGN/C0M8/R0EA,R0EB,R0EC  , ROEC 2 ,C PO ,C PI ,CP2,CP3, 
lCF0»CFl,CF2,CGDl,CG02tCGBl,CGB2tCGFl,CGF2 
IF     (IQ-1)     100,101,102 

100  IP=2 

GO  TO  200 

101  IP=1 

GO  TO  200 

102  1=1 
IP=2 

105  DO  195  J=IP,LL 

IF(LCUTU)}     195,110,195 
110    DRX=RX( J)-RX( I ) 

IF(DRX)     113,117,117 
113     IF(DRX+ROE)    195,195,120 
117    IFLDRX-ROE)     120,195,195 
120    DRY=RY( J )-RY( I ) 

IF(DRY)     123,127,127 
123    IF(DRY+ROE)     195,195,130 
127     IF(DRY-ROE)     130,195,195 
130    DRZ=RZ( J)-RZ(I ) 

IF(DRZ)  133,137,137 
133  IF(DRZ+ROE)  195,195,140 
137  IF(DRZ-ROE)  140,195,195 
140  DIST=DKX*DRX+DRY*DRY+DRZ*DRZ 

IF(DIST-R0E2)  150,195,195 
150  DIST=SQRT(DIST) 
160  IF(DIST-ROEK)  162,162,165 
162  FORCE=EXPlPFXA+PEXB*DIST ) 

GO    TO    180 
165    DFF=ROE-DIST 

IF(DFF-l.OE-lO)     195,195,167 
167    FORCE=(EXP(PAC+PEXB*DI ST)-PFPTC)/DFF 
180     IF(FM-FORCE)     190,190,195 
190    FOD=FORCE/DIST 

FA=FOD*DRX 

FX(J )  =  FX(J  )+FA 

FX( I )=FX( I )-FA 

FA=FOD*DRY 

FY(J)=FY(J)+FA 
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FY(I )=FY(I)-FA 
FA=FOD*DRZ 
FZ(J)=FZ(J)+FA 
FZ(I )=FZ( I J-FA 

CONTINUE 


195 

200  DO  300  1=1 
IF(LCUTU) 
205  IP=I+1 

DO  295  J=I 

IF(LCUT( J) 

210  DRX=RX(J)- 

IF(DRX)  21 

213  IF(DRX+ROE 

217  IF(DRX-ROE 

220  DRY=RY(J)- 

IF(DRY)  22 

223  IF(DRY+ROE 

227  IF(DRY-ROE 

230  DRZ=RZ(J)- 

IF(DRZ)  23 

233  IFIDRZ+RCE 

237  IF(DRZ-ROE 

240  DIST=DRX*D 

IF(DIST-RO 

250  DIST=SQRT( 

IF(DIST-RO 

255  IF(DIST-P.O 

260  FORCE=EXP( 

GO  TO  280 
265  FORCE=DIST 

GO  TO  280 
270  FORCE=EXP( 
280  IF(ABS(FOR 
FOD  =  FQRC 
FA=FOD*DRX 
FX( J)=FX( J 
FX{ I )=FX(1 
FA=FOD*DRY 
FY( J)=FY{ J 
FY(I J=FY(I 
FA=FOD*DRZ 
FZ( J)=FZ( J 
FZ( I )  =  FZ( I 
295  CONTINUE 
300  CONTINUE 
RETURN 
END 


P,LD 

)  300,205,300 


P,LL 
)  295, 
RX(I  ) 

3,217, 
C)  295 
C)  220 
RY(I) 

3,227, 
C)  295 

C)  230 
RZ(I  ) 

3,237, 
C)  295 
C)  240 
RX+DPY 
EC2)  2 
DIST) 

EA)  26 

EB)  26 
FXA+EX 


210,295 

217 

,295,220 
,295,295 

227 

,295,230 
,295,295 

237 

,295,240 

,295,295 

*DRY+DRZ*DRZ 

50,295,295 

0,255,255 
5,270,270 
B-DIST) 


*(DIST*CF2+CF1 )+CFO 

CGF1+CGB1*DIST )-EXP ( CGF2+CGB2*DI ST ) 

CE) .LE.FM)  GO  TO  295 

E/DIST 

)+FA 
)-FA 

)+FA 

)-FA 

)+FA 

)-FA 


SUBROUTINE  ENERGY 
C 

THIS  SUBROUTINE  CALCULATES  THE  MUTUAL  POTENTIAL  ENERGIES. 
THE  FIRST  HALF  DOES  THE  DYNAMICS  FOR  ATOM  #1 ;  THE  SECOND 
HALF  FOR  ALL  OTHERS. 
C 

COMMON/COM 1/RX( 1000) , RY( 1000 ), RZ (1000 ) ,LCUT(1000) , 
1LL,LD, ITYPE,NVAC 

C0MMCN/CQM5/RCE,R0E2,R0EM, EX A, EXB, PEXA , PEXB , FXA, PFX A, 
1IQ,TSAVE,BSAVE 

COMMON /C0M7/P  PTC,  T  POT,  PPE(  1000)  ,  TLP  E  ,  R  OEL  ,  ROcL  2,  ,MEW 

COMMON/COM 8/ROE A, ROE B,ROEC ,ROEC2,CPO,CP1,CP2,CP3, 
1CF0,CF1,CF2,CGC1,CGD2,CG31 ,CG32 ,CGF 1 ,CGF2 

IF    (IQ-1)     100,101,102 

IP=2 


100 
101 
102 
105 


GO  TO  200 

IP  =  1 

GO  TO  200 

1  =  1 

IP=2 

DO  595  J=IP,LL 

IF(LCUT(J)  )  595, 510,595 


61 


510 

513 
517 
520 

523 
527 
530 

533 
537 
540 

550 
560 
580 


595 
600 

200 

205 

210 

213 
217 
220 

223 
227 

230 


237 
240 

250 

255 
260 

265 

270 
280 


295 
300 


DRX=RX{ J)- 
I  F  (  D  R  X )  51 
IFCDRX+ROE 

IF(DRX-ROE 

DRY=RY( J)- 

IF(DRY)  52 

IF(DRY+ROE 

IF(DRY-ROE 

DRZ^=RZ(  J)- 

IF(DRZ)  53 

IFIDRZ+ROE 

IF(DRZ-RGE 

DIST=ORX*D 

IFIDIST-RO 

DIST=SQRT( 

POT=EXP(PE 

TPOT=TPOT+ 

PPE( 1)=PPE 

PPE{ J) =PPE 

CONTINUE 

CONTINUE 


RX(  I  ) 
3,517, 
)  595, 
)  520, 
RY(I  ) 
3,527, 
)  595, 
)  530, 
RZ(I  } 
3,537, 
)  595, 
)  540, 
RX+ORY 
E2)  55 
DIST) 
XA+PEX 
POT 

(I )+BS 
UJ  +  TS 


517 

595,520 
595,595 

527 

595,530 
595,595 

537 

595,540 
595, 595 
*DRY+DRZ*DRZ 
0,595,595 

B*DIST)-PPTC 

AVE*POT 
AVE*POT 


DO  3 
IF(L 
I  P  =  I 

DO  2 
IF(L 

DRX  = 
IF(D 
IFID 
IF(D 
DRY= 
IF(D 
IF(D 
IF(D 
DRZ  = 
IF(D 
IF  (D 
IF(D 
DIST 
IF(D 
DIST 
IF(D 
IF(D 
POT  = 
GO  T 
POT  = 
GO  -T 
POT  = 
TPOT 
SAVE 
PPE( 
PPE( 
CONT 
CO  NT 
RETU 
END 


00 
CU 
+  1 
95 
CU 
RX 
RX 
RX 
RX 
RY 
RY 
RY 
RY 
RZ 
RZ 
RZ 
RZ 
=  D 
IS 

=s 

IS 

IS 

EX 

0 

DI 

0 

EX 

=  T 

=  0 

I) 

J) 

IN 

IN 

RN 


I=IP,LD 
T(  I  )  )  300,205,300 


J  = 

T{J 
(J) 
)  2 

+  R0 
-RC 
(J) 
)  2 
+  R0 
-RO 
(J) 
)  2 
+  R0 
-RO 
RX* 
T-R 
QRT 
T-R 
T-R 
P(E 
280 
ST* 
280 
P(C 
POT 
.5* 
=  PP 
=  PP 
UE 
UE 


IP, LI 

))  295,210,295 

-RXU  ) 

13,217, 

EC)  295 

EC)  220,295,295 

-RY( I ) 

23 , 227 , 

EC)  295 

EC)  230 

-RZ(I ) 

33,237, 

EC)  295 


217 
,295,220 


227 

,295,230 

,295,295 


EC)  240 
DRX+DRY 
CEC2)  2 
(DIST) 
OEA)  26 
0E8)  26 
XA+EXB* 


,295,240 
,295,295 
^DRY+DRZ*DRZ 
50,295,295 

0,255,255 
5,270,270 
DIST) 


(DIST*(DIST*CP3+CP2)+CP1)+CP0 

GD1+CG31*DIST)-EXP(CGD2+CGB2*DIST) 

+  PCT 

POT 

E(l )+SAVE 

E(J)+SAVE 


C 

THIS 
TION 
C 


SUBROUTINE  PRINT 


SUBROUTINE 
AT  THE  TOP 


PRINTS  THE  HEADING  OF  ALL  PERTINENT  INFCRMA- 
OF  EACH  TIMESTEP  PRINTOUT. 


COMMON/COM1/RXU000)  ,RY(1000),RZ(1000)  ,LCUT(  1000)  , 
1LL,LD, ITYPEtNVAC 

C0MM0M/C0M2/IHK20),  I  H2  (  8  )  ,  IHS  (  10)  ,  I  HB  (  6  )  ,  I  HT  (  6)  , 
1TARGET(4)  ,TMAS,  BULLET  (4)  ,3  MAS-,  PLANE,  TEMP,  THERM 

COMMON /C0M3/KX I (  1000 ) , RY I {  1000 ) , RZ I ( 1000 )  , C VR  ,  E VR , 
1NT,TIME,DT,DTI ,  I  LAY 

COMMON /C0M4/IX, I Y,IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z, 
1IVACX,  1VACY, IVACZ 
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C0MM0N/C0M5/R0E,R0F2,R0EM, EX A, EX B, P EXA , P EXB , FXA, PFXA , 

HQ,TSAVE,eSAVE 
C0MM0N/C0M8/RQEA,R0EB,R0EC  ,R0EC2  ,CPO,CP  1  ,CP2  ,CP3  , 

1CF0,CF1,CF2,CGD1,CGD2,CGB1 , CGB2 , CGF1 ,  CG^2 
9710    FORMAT !40X,10A4, / , 28 X, 2  0A4 » / ) 
9720    F0RMAT(9H    TARGET    - ,4A4, 10HPRIMAR Y   -     ,4A4  ,  IX,  14HLATTI  CE 

1    UNIT    =,F7.4,4H    ANG) 
9730    FORMAT (4X,6HMASS    =, F7. 2, 13X, 6HMASS    = , F7. 2 , 9X , 14HLATT I C 

IE  TEMP  =F5.2,7H  DEG  K,,18H   THERMAL  CUTOFF  =,F5.2,3H  E 

1V/) 

9740  FORMAT (2H    (,A4,8H)     PLANE, ,18H       PRIMARY    ENERGY    =  , 

1  F5.2,21HKEV,   CRYSTAL  SIZE  (  ,12, 3H  X  ,12, 3H  X  ,I2,3H 
1  ),,  4X,  16HVACANCY  IN  SITE  ,  14/) 

9741  FORMAT (2H  (,A4,8H)  PLANE, , 1 8H   PRIMARY  ENERGY  =, 

1    F5.3,21HKCV,       CRYSYAL    SIZE     (     ,I2,3H    X    ,I2,3H    X    ,12,3 
1H    ),,    4X, 'INTERSTITIAL     ( ■ ,2 ( F5. 2, ' , ' ) , F5.2 , •     )     FROM 
1SITE     ■  ,14/) 

9742  F0RMAT(2H  (,A4,8H)  PLANE, ,18H   PRIMARY  ENERGY  =, 

1    F5.2,21HKEV,       CRYSTAL    SIZE     (     ,12, 3H    X     ,12, 3H    X     ,12, 3H 
1     ),,    4X,    20HREPLACEMENT     IN    SITE     ,     14/) 
9750    F0RMAT(30H    PRIMARY    START    POINT     ( LU )       X    =,F5.2,5H,     Y    =, 
1F5.2,5H,     Z    =,F5.2, 5X, 13, ■     LAYERS    ARE    FREE    TO    MOVE', 
110X,4HIQ    =,12/) 
9760    F0RMATC12H    POTENTIAL       ,6A4,3X,5HPEXA=, F9.5 , 2X, 5HPEXB=, 

1F9.5,2X,5HPFXA=.F9.5) 
9765    FORMAT ( 12X ,6 A4 , 3X, 5HEX A    =, F9 . 5 , 2X, 5HEXB    = , F9. 5, 2X, 5HFX 

1A    =,F9.5/) 
9770    FORMAT  ('     WHEN«,F8.4,  •     <    R    0,F8.4,  •     THE    MATCHING    POTEN 
1TIAL    PARAMETERS    ARE*,//,'  CPO    =',F10.3,',     CP1    =  • 

1F10.3,1,    CP2    =«,F10.3,»,    CP3    =«,F10.3,/,«  CFO    =  ■ 

1E10.3,',     CF1    =',E10.3,',     CF2    =',E10.3,//) 
9780    FORMAT! •     CUT-OFF    AT',F5.2,',    WHEN    R    >     ',F6.3,«    LU,     MOR 
1SE    POTENTIAL    PARAMETERS    ARE',     8A4,//,10X,'  CGD1     =', 

1F8.4,1,     CG02    =,,F8.4,1,     CGBL    =«,F8.4,»,    CGB2    =',F8.4, 
l'i    CGF1    =«,F8.4,»,    CGF2    =«,F3.4,//) 
9790    FORMAT!  10H    TIMESTEP     , 1 4 , 22X , 6HDTI     =     ,     F5.4,     5H    LU, 
1,22H       ELAPSED    TIME     (SEC)     =,     E10.4, •,     NEXT    TIMESTEP     IS 
1=' ,E10.4/) 

WRITE     (    6,9710)     IHS,IH1 
WRITE     (    6,9720)     T ARGET , BU LLET , CVR 
WRITE     (     6,9730}     TMAS , BMAS , TEMP , TH ERM 
GO    TO     (401,402,403,402),     ITYPE 

401  WRITE     (     6,9740)     PLANE , E VR ,  IX, I Y , I Z , NV AC 
GO    TO    405 

402  WRITE     (    6,9741)     PLANE , E VR , IX , I Y , I Z , D1X, Dl Y, D 1Z, NV AC 
GO    TO    405 

403  WRITE     (     6,9742)     P LANE ,E VR , I  X , I Y , I Z , NVAC 

405       WRITE     (    6,9750)     RXI ! 1)  , RY  I  ( 1 )  ,  R Z I (  1 )  ,  IL AY ,  I  Q 
WRITE     (     6,9760)     I HB  ,  PEXA , P EXB , PFX A 
WRITE     !     6,9765)     IHT , EXA ,E XB ,FXA 

WRITE     !    6,9770)    ROE  A, ROEB , CPO , CP1 , CP2, CP3 , C FO,CF 1 , CF2 
WRITE     I    6,9780)     ROEC , ROEB , I  HZ , CGD1 , CGD2 , CGB1 , CGB2 , 
1CGF1,CGF2 

WRITE     !    6,9790)     NT , DT I , T I  ME, DT 
RETURN 
END 

BLOCK    DATA 
DIMENSIONING    OF    VARIABLES    USED    IN    COMMON 

COMMON/COMl/RXtlOOOJ ,RY( 1000 ), RZ ( 10  00 ) ,LCUT( 1000) , 
1LL,LD, ITYPE, NVAC 
DATA    R X/ 1 000*0. 0 /, R Y/ 1000* 0. 0/,RZ/l 000*0. 0/,LCUT/ 1000* 
CCMM0N/C0M3/RXI ( 1000) , RY I ( 1000),RZI( 10  00 ) ,  C VR , E VR , 
INT, TIME, DT, DTI , ILAY 
DATA    RXI / 1 000*0. 0/, RYI/1 000*0. 0/,RZI/ 1000*0.0/ 
C0MM0N/C0M6/FX ( 1 000 )  , FY ( 1000  ) , FZ ( 1 000 ) , PAC , P FDTC , FM 
DATA    F X/ 1000*0. 0 /, F Y/ 1 000* 0.0/, FZ/ 1000* 0.0/ 
END 
//GO.FT06F001    DD    SPACE= ( CYL , ! 1 0 , 2 ) , RLS E) 
//GO.SYSIM    DD    * 
CRYSTAL-1968    MODIFIED    TO    DEAL    WITH    VACANCIES    AND    INTERS.TITI 
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(  GIRIFALCO — WEIZER  POTENTIAL  )   .99060  1.41160  3.03200  3.4 

ARGON  39.948     9.33    -5.60  ARGON-TUNGSTEN 

TUNGSTEN        183.86    11.30    -7.50  TUNGSTEN-TUNGSTEN 

BODY  CENTERED  CUBIC,  (100)  ORIENTATION  100     10  10  10 

100        5  4   51  0.7  -0.7   0.0  5  1 
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